This script fits state-space models to TRPD grassland bird survey data to assess species-level population trends.
knitr::opts_chunk$set(echo = TRUE)
library(here)
library(png)
library(jagsUI)
library(ggplot2)
library(ggthemes)
library(tidyverse)
library(kableExtra)
source(here("Scripts", "1_Data_Import.R"))
source(here("Scripts", "2_Spatial_Data_Import.R"))
source(here("Scripts", "3_Data_Cleaning.R"))
source(here("Scripts", "4_Covariate_Data_Prep.R"))
source(here("Scripts", "5_Prep_Species.R"))# Set random seed
set.seed(42563)
# Identify species and transsects for analysis.
species_to_run <- prairie_sp
transects_to_sum <- c("CHOA01", "CHOA02", "CHOA09", "CHOA14")
# ID tag for this run and create associated output directory
tag <- "20231128_plus1_to_counts"
mainDir <- here("Results", "991_bird_trends_species")
dir.create(file.path(mainDir, tag), showWarnings = TRUE)
outputDir <- here("Results", "991_bird_trends_species", tag)
subdirectories <- c("annual_n_estimates", "annual_n_plots", "diagnostic_plots", "species_level_estimates", "model_outputs", "synthesized_results", "Rhats", "model_parameters")
lapply(file.path(mainDir, tag, subdirectories), function(x) if (!dir.exists(x)) dir.create(x))
adjust_counts_by_1 <- TRUE
write.table(transects_to_sum, file = here("Results", "991_bird_trends_species", tag, "model_parameters", "transects_summed.txt"))transect_birds <- bird.tblBirds_Transects
# Merge BB and OA surveys that are actually the same transect
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB02")] <- "CHOA02"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB09")] <- "CHOA09"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB14")] <- "CHOA14"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB16")] <- "CHOA14" # shouldn't be any...
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB01")] <- "CHOA01"
# Filter bird data to only specified transects
transect_birds <- transect_birds %>%
left_join(bird.tblSurveyDefn, by = "SurveyCode") %>%
filter(SurveyCode %in% transects_to_sum) %>%
mutate(
survey_string = paste0(SurveyCode, "-", SurveyDate), # unique code for survey based on area, site #, and date
year = sub("-.*", "", SurveyDate), # convert survey date to just year
year = as.numeric(year)
) # year to numeric
# Identify days and years in which all transects in the list were actually surveyed
survey_dates <- transect_birds %>%
group_by(SurveyDate, SurveyCode) %>%
summarize(survey_count = length(unique(SurveyDate))) %>%
pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
ungroup() %>%
filter(complete.cases(.)) %>%
select(SurveyDate)
survey_complete_years <- transect_birds %>%
group_by(year, SurveyCode) %>%
summarize(survey_count = length(unique(year))) %>%
pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
ungroup() %>%
filter(complete.cases(.)) %>%
select(year)
# Re-filter bird data to only dates with complete surveys across all transects. This version forces all transects to have been run on the same day.
transect_birds <- transect_birds %>%
filter(SurveyDate %in% survey_dates$SurveyDate)
# In years surveyed more than once, pick just the last survey from that year.
surveys_all <- transect_birds %>%
group_by(year, SurveyCode, survey_string) %>%
summarize() %>%
mutate(surveycode_year = paste0(SurveyCode, "-", year))
surveys_selected <- surveys_all[!duplicated(surveys_all$surveycode_year, fromLast = TRUE), ] %>%
ungroup() %>%
select(survey_string)
# Re-filter bird data based on previous step to only include one survey day each year.
transect_birds <- transect_birds %>%
filter(survey_string %in% surveys_selected$survey_string)
# Only include counts > 0 (remove rows where bird only detected outside survey area).
transect_birds <- transect_birds %>%
filter(BirdCountIn > 0)
# Get species list from data and identify species observed at least 5 times.
sp_detected <- unique(transect_birds$BirdCode)
sp_detections <- transect_birds %>%
group_by(BirdCode) %>%
filter(BirdCode %in% prairie_sp) %>%
summarise(
count_indiv = sum(BirdCountIn),
n_surveys = n(),
n_yrs = length(unique(year))
)
sp_detected_above_threshold <- sp_detections %>%
filter(count_indiv >= 5) %>%
pull(BirdCode)
# Re-filter species to run so that it only includes species actually detected in selected surveys
species_to_run <- intersect(species_to_run, sp_detected)
# Also create version with only species detected more than 5 times.
species_to_run_above_threshold <- intersect(species_to_run, sp_detected_above_threshold)
# Save these species list for use in subsequent analyses
saveRDS(sp_detected_above_threshold, here("Results", "991_bird_trends_species", tag, "sp_detected_above_threshold.RDS"))
saveRDS(species_to_run_above_threshold, here("Results", "991_bird_trends_species", tag, "species_to_run_above_threshold.RDS"))This code creates a function for preparing time series data for species of interest in transects of interest and then fits a state-space model (assuming density independence) to the data.
# Function to prepare time series data for species and transects of interest
prep_and_fit_CHBB <- function(species_of_interest = species_to_run, transects_to_assess = transects_to_sum) {
transect_birds <- bird.tblBirds_Transects
# Merge BB and OA surveys that are actually the same transect
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB02")] <- "CHOA02"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB09")] <- "CHOA09"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB14")] <- "CHOA14"
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB16")] <- "CHOA14" # shouldn't be any...
transect_birds$SurveyCode[which(transect_birds$SurveyCode == "CHBB01")] <- "CHOA01"
# Filter bird data to only those transects
transect_birds <- transect_birds %>%
left_join(bird.tblSurveyDefn, by = "SurveyCode") %>%
filter(SurveyCode %in% transects_to_assess) %>%
mutate(
survey_string = paste0(SurveyCode, "-", SurveyDate), # Unique code for survey based on area, site #, and date
year = sub("-.*", "", SurveyDate), # Convert survey date to just year
year = as.numeric(year)
)
# Identify days and years in which all transects in the list were actually surveyed
survey_dates <- transect_birds %>%
group_by(SurveyDate, SurveyCode) %>%
summarize(survey_count = length(unique(SurveyDate))) %>%
pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
ungroup() %>%
filter(complete.cases(.)) %>%
select(SurveyDate)
survey_complete_years <- transect_birds %>%
group_by(year, SurveyCode) %>%
summarize(survey_count = length(unique(year))) %>%
pivot_wider(names_from = SurveyCode, values_from = survey_count) %>%
ungroup() %>%
filter(complete.cases(.)) %>%
select(year)
# Re-filter bird data to only dates with complete surveys across all transects.This version forces all transects to have been run on the same day.
transect_birds <- transect_birds %>%
filter(SurveyDate %in% survey_dates$SurveyDate)
# In years surveyed more than once, pick just the last survey from that year
surveys_all <- transect_birds %>%
group_by(year, SurveyCode, survey_string) %>%
summarize() %>%
mutate(surveycode_year = paste0(SurveyCode, "-", year))
surveys_selected <- surveys_all[!duplicated(surveys_all$surveycode_year, fromLast = TRUE), ] %>%
ungroup() %>%
select(survey_string)
# Re-filter bird data based on the previous step to only include one survey day each year
transect_birds <- transect_birds %>%
filter(survey_string %in% surveys_selected$survey_string)
# Create df of years with complete data
survey_years <- unique(transect_birds$year) %>%
sort() %>%
as.data.frame()
survey_years <- rename(survey_years, year = .)
# Count the number of species of interest across all transects
sp_counts <- transect_birds %>%
filter(BirdCode == species_of_interest) %>%
group_by(year) %>%
summarize(total_count = sum(BirdCountIn))
# Zero fill for surveyed years
sp_counts <- sp_counts %>%
right_join(survey_years) %>%
replace_na(list(total_count = 0)) %>%
arrange(year)
# This version of state-space model can't be fit when there are zeros in the data (populations that fit zero have no means of recovering/growing). Adding 1 to everything for now...
sp_counts <- sp_counts %>%
mutate(
total_count_orig = total_count
)
if (adjust_counts_by_1 == TRUE) {
sp_counts <- sp_counts %>%
mutate(total_count = total_count + 1)
}
# Identify information about years for model
yr_min <- min(sp_counts$year)
yr_max <- max(sp_counts$year)
year_range <- seq(yr_min, yr_max, 1) # Create a vector of survey years
n_yrs <- length(year_range) # Counter for how many years of data
year_range_df <- year_range %>%
as.data.frame()
year_range_df <- rename(year_range_df, year = `.`)
# Add NAs for un-surveyed years
sp_counts <- sp_counts %>%
right_join(year_range_df) %>%
arrange(year)
#--------------------------
# JAGS code for Bayesian hierarchical exponential population growth state space model
{
sink("Scripts/JAGS_code/991_bird_trends_species_version_for_MOU.jags")
cat("
model {
# Specify prior values
init_lnN ~ dnorm(ln1, 1)
lnN[1] <- init_lnN
r_0 ~ dnorm(0,1) # mean intrinsic growth rate
sigma_r ~ dunif(0, 20) # annual SD in true pop growth rate
obs_sd ~ dunif(0,20) # annual SD in observation (measurement) error
# convert SD to precision (tau=1/SD^2)
sigma_tau <- pow(sigma_r, -2)
obs_tau <- pow(obs_sd, -2)
# Likelihood
for (t in 1:(n_yrs-1)){
lnN[t+1] ~ dnorm(lnN[t] + r_0, sigma_tau)
}
# Observation process
for (t in 1:n_yrs) {
lnCount[t] ~ dnorm(lnN[t], obs_tau)
N_est[t] <- exp(lnN[t]) # convert log(N) back to real scale
} # close t loop
# Derived parameters
# Geometric mean rates of change (%/year between first and last year)
ch <- N_est[n_yrs] / N_est[1]
trend <- (ch^(1/(n_yrs-1)))-1
# Trend from lognormal slope
slope_trend <- ( n_yrs * sum(1:n_yrs * lnN) - sum(1:n_yrs) * sum(lnN) ) / (n_yrs * sum((1:n_yrs)^2) - sum(1:n_yrs)^2)
} # End jags model
", fill = TRUE)
sink()
}
# Bundle data, identify parameters to save, submit to jags
TRPD_data <- list(
n_yrs = n_yrs, lnCount = log(sp_counts$total_count),
ln1 = log(sp_counts$total_count[1])
)
parameters <- c("r_0", "sigma_r", "obs_sd", "N_est", "trend", "slope_trend")
TRPD1 <- jagsUI(TRPD_data,
inits = NULL, parameters, "Scripts/JAGS_code/991_bird_trends_species_version_for_MOU.jags",
n.adapt = 1000, n.chains = 3, n.thin = 10,
n.iter = 300000, n.burnin = 50000, parallel = TRUE
)
#---------------------------------------------
# Save model output. Models are too big to save everything, not including MCMC chains.
mod_to_save <- TRPD1[which(sapply(TRPD1, class) != "mcmc.list")]
saveRDS(mod_to_save, file = paste0(outputDir, "/model_outputs/", species_of_interest, "_model_outputs.RDS"))
# Assess convergence by verifying all R-hat values < 1.1
max.rhat <- function(x) {
max(max_rhat <- c(sapply(x$Rhat, max, na.rm = TRUE)))
}
max.rhat.run <- max.rhat(TRPD1) # r.hat <1.1 suggests good convergence
saveRDS(TRPD1$Rhat, file = paste0(outputDir, "/Rhats/", species_of_interest, "_rhats.RDS"))
# Summarize r estimates
r_est <- paste0("r: ", round(as.numeric(TRPD1$mean["r_0"]), 2), " (95% CI: ", round(as.numeric(TRPD1$q2.5["r_0"]), 2), "-", round(as.numeric(TRPD1$q97.5["r_0"]), 2), ")")
# Summarize percent of r_0 posterior distribution above 0
f_r0 <- round(as.numeric(TRPD1$f["r_0"]), 2)
# Data frame with model estimates
sp_values <- data.frame(
sp = species_of_interest,
r0_mean = TRPD1$mean$r_0,
r0_sd = TRPD1$sd$r_0,
r0_median = TRPD1$q50$r_0,
r0_LCI_95 = TRPD1$q2.5$r_0,
r0_UCI_95 = TRPD1$q97.5$r_0,
r0_LCI_50 = TRPD1$q25$r_0,
r0_UCI_50 = TRPD1$q75$r_0,
f_r0 = round(as.numeric(TRPD1$f["r_0"]), 2),
sigma_r_mean = TRPD1$mean$sigma_r,
sigma_r_sd = TRPD1$sd$sigma_r,
sigma_r_median = TRPD1$q50$sigma_r,
sigma_r_LCI_95 = TRPD1$q2.5$sigma_r,
sigma_r_UCI_95 = TRPD1$q97.5$sigma_r,
sigma_r_LCI_50 = TRPD1$q25$sigma_r,
sigma_r_UCI_50 = TRPD1$q75$sigma_r,
f_sigma_r = round(as.numeric(TRPD1$f["sigma_r"]), 2),
obs_sd_mean = TRPD1$mean$obs_sd,
obs_sd_sd = TRPD1$sd$obs_sd,
obs_sd_median = TRPD1$q50$obs_sd,
obs_sd_LCI_95 = TRPD1$q2.5$obs_sd,
obs_sd_UCI_95 = TRPD1$q97.5$obs_sd,
obs_sd_LCI_50 = TRPD1$q25$obs_sd,
obs_sd_UCI_50 = TRPD1$q75$obs_sd,
f_obs_sd = round(as.numeric(TRPD1$f["obs_sd"]), 2),
max.rhat.run = max.rhat(TRPD1),
method = "SSMr"
)
write.csv(sp_values, file = paste0(outputDir, "/species_level_estimates/", species_of_interest, "_species_level_estimates.csv"))
# Add estimates to simple dataframe
sp_counts_w_ests <- sp_counts
sp_counts_w_ests <- sp_counts_w_ests %>%
mutate(
N_est_mean = TRPD1$mean$N_est,
N_est_sd = TRPD1$sd$N_est,
N_est_median = TRPD1$q50$N_est,
N_est_LCI_95 = TRPD1$q2.5$N_est,
N_est_UCI_95 = TRPD1$q97.5$N_est,
N_est_LCI_50 = TRPD1$q25$N_est,
N_est_UCI_50 = TRPD1$q75$N_est,
sp = species_of_interest,
method = "SSMr"
)
write.csv(sp_counts_w_ests, file = paste0(outputDir, "/annual_n_estimates/", species_of_interest, "_annual_N_ests.csv"))
# Plot posterior draws
output_path <- here("Results", "991_bird_trends_species", tag, "diagnostic_plots", paste0(species_of_interest, "_diagnostic_plot.png"))
png(output_path)
par(mfrow = c(2, 1))
plot(TRPD1$sims.list$obs_sd, TRPD1$sims.list$sigma_r,
ylab = "Process error",
xlab = "Observation error", cex = 0.3
)
plot(TRPD1$sims.list$r_0, TRPD1$sims.list$sigma_r,
ylab = "Process error",
xlab = "Process mean", cex = 0.3
)
dev.off()
if (adjust_counts_by_1 == TRUE) {
sp_counts_w_ests2 <- sp_counts_w_ests %>%
mutate(
N_est_LCI_95 = N_est_LCI_95 - 1,
N_est_UCI_95 = N_est_UCI_95 - 1,
N_est_mean = N_est_mean - 1,
total_count = total_count - 1
) %>%
mutate(
N_est_LCI_95 = ifelse(N_est_LCI_95 < 0, 0, N_est_LCI_95),
N_est_UCI_95 = ifelse(N_est_UCI_95 < 0, 0, N_est_UCI_95),
N_est_mean = ifelse(N_est_mean < 0, 0, N_est_mean)
)
} else {
sp_counts_w_ests2 <- sp_counts_w_ests
}
# Make time series plot
plot <- ggplot(sp_counts_w_ests2, aes(x = year)) +
geom_ribbon(aes(ymin = N_est_LCI_95, ymax = N_est_UCI_95), fill = "grey90", alpha = .6) +
geom_line(aes(y = N_est_mean), color = "black", alpha = .7) +
geom_point(aes(y = N_est_mean), color = "black", alpha = .7) +
geom_point(aes(y = total_count), color = "maroon", shape = 1) +
annotate(geom = "text", label = r_est, x = min(sp_counts_w_ests$year) + 5, y = max(sp_counts_w_ests$total_count), size = 3) +
scale_x_continuous(breaks = seq(1960, 2021, by = 10)) +
theme_clean() +
labs(
title = paste0(bird.tblBirdSpecies$BirdName[which(bird.tblBirdSpecies$BirdCode == species_of_interest)], " population dynamics"),
subtitle = paste0("TRPD transects- ", paste(transects_to_assess, sep = "", collapse = ", "), "\nModeled estimate with 95% CI and raw counts \n", r_est),
x = "Year",
y = "Total count across transects"
) +
theme(
plot.title = element_text(size = 11, face = "bold"),
plot.subtitle = element_text(size = 10, face = "italic"),
legend.title = element_text(size = 10, face = "plain")
)
ggsave(plot = plot, filename = paste0(species_of_interest, "_plot.png"), device = "png", path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), width = 5, height = 3)
}# # fit state space models for species of interest
# results <- lapply(species_to_run, prep_and_fit_CHBB)
# names(results) <- species_to_run
species_still_to_run <- list.files(
path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), pattern = NULL, all.files = FALSE,
full.names = FALSE, recursive = FALSE,
ignore.case = FALSE, include.dirs = FALSE, no.. = FALSE
)
species_still_to_run <- gsub("\\_plot.png", "", species_still_to_run)
species_still_to_run <- unique(species_still_to_run)
species_still_to_run <- setdiff(species_to_run, species_still_to_run)
library(parallel)
numCores <- detectCores() # get the number of cores available
if (length(species_still_to_run) > 0) {
lapply(species_still_to_run, prep_and_fit_CHBB)
# mclapply(species_to_run, prep_and_fit_CHBB, mc.cores = numCores)
}Load model results
species_level_estimates <- list.files(here("Results", "991_bird_trends_species", tag, "species_level_estimates"), pattern = "species_level_estimates.csv", full.names = TRUE) %>%
lapply(read.csv) %>%
bind_rows() %>%
select(-X) %>%
filter(sp %in% species_to_run_above_threshold) %>%
mutate(
r0_mean = round(r0_mean * 100, 1),
r0_LCI_95 = round(r0_LCI_95 * 100, 1),
r0_UCI_95 = round(r0_UCI_95 * 100, 1)
)The models allow us to estimate the contribution of both observer error (undercounting/overcounting) and process error (real variability in growth rates) to the observed annual variability in populations.
error_parsing <- ggplot(data = species_level_estimates, aes(x = sigma_r_mean, y = obs_sd_mean, label = sp)) +
geom_linerange(aes(xmin = sigma_r_LCI_50, xmax = sigma_r_UCI_50)) +
geom_linerange(aes(ymin = obs_sd_LCI_50, ymax = obs_sd_UCI_50)) +
geom_label(size = 2.4, label.padding = unit(.075, "lines")) +
theme_clean() +
labs(
x = "Process error (Growth variate variability)",
y = "Observer error"
)
error_parsing
error_parsing_repel <- ggplot(data = species_level_estimates, aes(x = sigma_r_mean, y = obs_sd_mean, label = sp)) +
ggrepel::geom_label_repel(size = 2.4, box.padding = 0.025, label.padding = .075, arrow = arrow(length = unit(0.01, "npc"), type = "open"), max.time = 20, max.iter = 100000) +
theme_clean() +
labs(
x = "Process error (Growth variate variability)",
y = "Observer error"
)
error_parsing_repel
ggsave(plot = error_parsing, filename = "error_parsing_plot.png", device = "png", path = here("Results", "991_bird_trends_species", tag, "synthesized_results"), width = 7, height = 7)
ggsave(plot = error_parsing_repel, filename = "error_parsing_repel_plot.png", device = "png", path = here("Results", "991_bird_trends_species", tag, "synthesized_results"), width = 7, height = 7)files <- list.files(path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), full.names = TRUE)
files_names <- list.files(path = here("Results", "991_bird_trends_species", tag, "annual_n_plots"), full.names = FALSE)
files_names <- gsub("\\_plot.png", "", files_names)
files <- files[which(files_names %in% species_to_run_above_threshold)]
files_names <- files_names[which(files_names %in% species_to_run_above_threshold)]
for (i in 1:length(files)) {
cat(paste0("#### ", bird.tblBirdSpecies.ebird$common_name[which(bird.tblBirdSpecies.ebird$BirdCode == files_names[i])], "\n"))
cat(paste0("\n\n"))
cat(paste0("**Annual trend (% change/year), 1978-2023:** ", species_level_estimates$r0_mean[which(species_level_estimates$sp == files_names[i])], "% (95% CI: ", species_level_estimates$r0_LCI_95[which(species_level_estimates$sp == files_names[i])], "% to ", species_level_estimates$r0_UCI_95[which(species_level_estimates$sp == files_names[i])], "%)", "\n\n"))
}Annual trend (% change/year), 1978-2023: 0.9% (95% CI: -1.5% to 3.3%)
Annual trend (% change/year), 1978-2023: 2% (95% CI: -2.6% to 7.4%)
Annual trend (% change/year), 1978-2023: -0.6% (95% CI: -5.1% to 4.6%)
Annual trend (% change/year), 1978-2023: -0.2% (95% CI: -20.7% to 20.5%)
Annual trend (% change/year), 1978-2023: -0.2% (95% CI: -11.6% to 11.5%)
Annual trend (% change/year), 1978-2023: 1.1% (95% CI: -7.6% to 10.8%)
Annual trend (% change/year), 1978-2023: 0.5% (95% CI: -1.9% to 3.3%)
Annual trend (% change/year), 1978-2023: -2.4% (95% CI: -10.6% to 5.1%)
Annual trend (% change/year), 1978-2023: 3.8% (95% CI: -7% to 15.5%)
Annual trend (% change/year), 1978-2023: 1.5% (95% CI: -1.3% to 5.1%)
Annual trend (% change/year), 1978-2023: 2.4% (95% CI: -3.6% to 9.1%)
Annual trend (% change/year), 1978-2023: 1% (95% CI: -4.8% to 6.7%)
Annual trend (% change/year), 1978-2023: 6.5% (95% CI: -0.8% to 13.9%)
Annual trend (% change/year), 1978-2023: 2.3% (95% CI: -9.2% to 13%)
Annual trend (% change/year), 1978-2023: 0.8% (95% CI: -4.7% to 6.9%)
Annual trend (% change/year), 1978-2023: 5.2% (95% CI: -6.8% to 16.9%)
Annual trend (% change/year), 1978-2023: 1% (95% CI: -3.6% to 5.8%)
Annual trend (% change/year), 1978-2023: -2.3% (95% CI: -10% to 5.6%)
Annual trend (% change/year), 1978-2023: 0.5% (95% CI: -5.4% to 6.2%)
Annual trend (% change/year), 1978-2023: 0.9% (95% CI: -1.5% to 3.9%)
Annual trend (% change/year), 1978-2023: 1.2% (95% CI: -11.2% to 13.9%)
Annual trend (% change/year), 1978-2023: -1.2% (95% CI: -7.7% to 5.1%)
Annual trend (% change/year), 1978-2023: 2.7% (95% CI: -5.5% to 11.6%)
Annual trend (% change/year), 1978-2023: -1.4% (95% CI: -13.9% to 10.9%)
Annual trend (% change/year), 1978-2023: 5.7% (95% CI: -3.9% to 15.4%)
Annual trend (% change/year), 1978-2023: -2.2% (95% CI: -9.1% to 6.5%)
Annual trend (% change/year), 1978-2023: -1.7% (95% CI: -9.6% to 5.5%)
Annual trend (% change/year), 1978-2023: 0.3% (95% CI: -4.6% to 4.8%)
Print session info:
## R version 4.2.2 (2022-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22621)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8
## [3] LC_MONETARY=English_United States.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_United States.utf8
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] jagsUI_1.5.2 png_0.1-8 transformr_0.1.3 gifski_1.12.0-2
## [5] gganimate_1.0.8 ggspatial_1.1.9 ggrepel_0.9.4 RVAideMemoire_0.9-83-7
## [9] pairwiseAdonis_0.4.1 cluster_2.1.4 goeveg_0.6.5 vegan_2.6-4
## [13] lattice_0.20-45 permute_0.9-7 mapview_2.11.2 ggtext_0.1.2
## [17] ratelimitr_0.4.1 rvest_1.0.3 trelliscopejs_0.2.6 plotly_4.10.3
## [21] auk_0.7.0 readxl_1.4.3 kableExtra_1.3.4 ggthemes_5.0.0
## [25] forcats_1.0.0 stringr_1.5.1 purrr_1.0.2 readr_2.1.4
## [29] tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.4 tidyverse_2.0.0
## [33] lubridate_1.9.3 dplyr_1.1.4 DT_0.31 sf_1.0-14
## [37] rgdal_1.6-4 sp_2.1-2 RODBC_1.3-23 here_1.0.1
##
## loaded via a namespace (and not attached):
## [1] spam_2.10-0 uuid_1.1-1 backports_1.4.1
## [4] Hmisc_5.1-1 systemfonts_1.0.5 plyr_1.8.9
## [7] lazyeval_0.2.2 splines_4.2.2 crosstalk_1.2.1
## [10] leaflet_2.2.1 unmarked_1.3.2 leafpop_0.1.0
## [13] digest_0.6.33 htmltools_0.5.7 leaflet.providers_2.0.0
## [16] fansi_1.0.5 magrittr_2.0.3 checkmate_2.3.0
## [19] tzdb_0.4.0 vroom_1.6.4 svglite_2.1.2
## [22] timechange_0.2.0 lpSolve_5.6.19 prettyunits_1.2.0
## [25] ggh4x_0.2.6 colorspace_2.1-0 leafem_0.2.3
## [28] textshaping_0.3.7 xfun_0.41 crayon_1.5.2
## [31] jsonlite_1.8.7 brew_1.0-8 glue_1.6.2
## [34] gtable_0.3.4 webshot_0.5.5 maps_3.4.1.1
## [37] scales_1.3.0 DBI_1.1.3 Rcpp_1.0.11
## [40] htmlTable_2.4.2 viridisLite_0.4.2 xtable_1.8-4
## [43] progress_1.2.3 gridtext_0.1.5 units_0.8-5
## [46] foreign_0.8-83 bit_4.0.5 proxy_0.4-27
## [49] mclust_6.0.1 dotCall64_1.1-1 Formula_1.2-5
## [52] stats4_4.2.2 htmlwidgets_1.6.4 epuRate_0.1
## [55] httr_1.4.7 wk_0.9.1 ellipsis_0.3.2
## [58] pkgconfig_2.0.3 farver_2.1.1 nnet_7.3-18
## [61] sass_0.4.7 utf8_1.2.4 tidyselect_1.2.0
## [64] labeling_0.4.3 rlang_1.1.2 later_1.3.1
## [67] munsell_0.5.0 cellranger_1.1.0 tools_4.2.2
## [70] cachem_1.0.8 cli_3.6.1 generics_0.1.3
## [73] evaluate_0.23 fastmap_1.1.1 yaml_2.3.7
## [76] ragg_1.2.6 knitr_1.45 bit64_4.0.5
## [79] s2_1.1.4 satellite_1.0.4 pbapply_1.7-2
## [82] nlme_3.1-160 mime_0.12 xml2_1.3.5
## [85] compiler_4.2.2 rstudioapi_0.15.0 e1071_1.7-13
## [88] tweenr_2.0.2 bslib_0.6.1 stringi_1.8.2
## [91] highr_0.10 fields_15.2 Matrix_1.6-4
## [94] classInt_0.4-10 commonmark_1.9.0 markdown_1.12
## [97] vctrs_0.6.5 pillar_1.9.0 lifecycle_1.0.4
## [100] jquerylib_0.1.4 MCMCvis_0.16.3 data.table_1.14.8
## [103] cowplot_1.1.1 raster_3.6-26 httpuv_1.6.12
## [106] R6_2.5.1 promises_1.2.1 gridExtra_2.3
## [109] KernSmooth_2.23-20 rjags_4-15 codetools_0.2-18
## [112] MASS_7.3-58.1 assertthat_0.2.1 rprojroot_2.0.4
## [115] withr_2.5.2 autocogs_0.1.4 mgcv_1.8-41
## [118] hms_1.1.3 terra_1.7-55 rpart_4.1.19
## [121] grid_4.2.2 rebird_1.3.0 coda_0.19-4
## [124] class_7.3-20 rmarkdown_2.25 DistributionUtils_0.6-1
## [127] shiny_1.8.0 base64enc_0.1-3
By Sam Safran